function [p_s, PM_s, ...
         rho_qm_ig_tau_ig, rho_qm_ig_taustar_ig, rho_qm_ig_a_ig, rho_qm_ig_zstar_ig, rho_qm_ig_astar_ig, rho_qm_ig_delta_ig, ...
         rho_pm_ig_tau_ig, rho_pm_ig_taustar_ig, rho_pm_ig_a_ig, rho_pm_ig_zstar_ig, rho_pm_ig_astar_ig, rho_pm_ig_delta_ig, ...
         rho_px_ig_tau_ig, rho_px_ig_taustar_ig, rho_px_ig_a_ig, rho_px_ig_zstar_ig, rho_px_ig_astar_ig, rho_px_ig_delta_ig, ...
         rho_pstar_ig_tau_ig, rho_pstar_ig_taustar_ig, rho_pstar_ig_a_ig, rho_pstar_ig_zstar_ig, rho_pstar_ig_astar_ig, rho_pstar_ig_delta_ig, ...
         rho_rm_ig_tau_ig, rho_rm_ig_taustar_ig, rho_rm_ig_a_ig, rho_rm_ig_zstar_ig, rho_rm_ig_astar_ig, rho_rm_ig_delta_ig] ...   
         = compute_equilibrium_foa_full( delta_ig, a_star_ig, z_star_ig, a_ig, aM_g, AM_s, AD_s, betaT, beta_s, alphaL_s, alphaI_s, alpha_ks, barL_rs, Z_rs, D, ... 
         tau_ig, tau_star_ig, m_ig_0, E_s_0, p_s_0, Dsg, DX, DM, omega_star, sigma, eta, kappa, sigma_star, nu, epsilon, gammaQ, gammaX, gammaM, rhoX, rhoM, tol, adj_inner, adj_outter, ...
         ind_Mshifts_vec, ind_Xshifts_vec, ind_Msample_ig_vec, ind_Xsample_ig_vec)

%This script computes the Jacobian matrix of the model predictions with respect to changes in tariffs
%Corresponding author: Rodrigo Adao
%Date: 09/11/2024

%% Preliminaries

D_s_last = E_s_0;
[R, S] = size(barL_rs);
[N, G] = size(a_ig);

%Inital product imports for correcting parameters
m_g0 = sum( m_ig_0, 1)'; 
%m_s0 = Dsg*m_g0;

%Compute parameters that do not change across iterations

    %terms for variety-level imports
    varphi_ig = ( a_ig.^(1+gammaM*omega_star) ).*( z_star_ig.^( gammaM*(1-sigma) ) ).*( (1 + tau_ig).^( gammaM*(1-sigma) ) );
    varphi_ig = varphi_ig.^( 1/(1+gammaM*omega_star*sigma) ) ;

    tilde_varphiR_ig = ( a_ig.^(1+omega_star) ).*( z_star_ig.^( 1 - gammaM*sigma ) ).*( (1 + tau_ig).^( - gammaM*sigma*(1+omega_star) ) );
    tilde_varphiR_ig = tilde_varphiR_ig.^( 1/(1+gammaM*omega_star*sigma) ) ;
    
    varphiR_ig = tau_ig.*tilde_varphiR_ig;

    %terms for product-level imports
    varphi_g = sum( varphi_ig , 1)';
    varphiR_g = sum( varphiR_ig , 1)';

    mu_g = ( aM_g.^( (1+gammaM*omega_star)/(1+gammaM*omega_star*eta) ) ).*( varphi_g.^( ( (1-eta)*(1+gammaM*omega_star*sigma) )/( (1-sigma)*(1+gammaM*omega_star*eta) ) ) ); 
    mu_gs_aux = ( ( Dsg*m_g0 )'*Dsg )'; %Wlog assign cost of one to products in sectors without impots 
    mu_g(mu_gs_aux<tol) = 1;
    mu_s = Dsg*mu_g;
    
    %terms for sector-level imports
    zeta_s = ( ( AM_s.^(gammaM*omega_star) ).*( mu_s.^( ( 1+gammaM*omega_star*eta )/(1-eta) ) ) ).^(1/(1+gammaM*omega_star*kappa)) ;
    
    %terms for production
    alphaK_s = 1 - alphaL_s - alphaI_s;
    alphaK_rs = repmat(alphaK_s',R,1);
    alphaI_rs = repmat(alphaI_s',R,1);
    alphaL_rs = repmat(alphaL_s',R,1);

    %terms for exports
    deltaX_ig_term = a_star_ig.*( (1 + tau_star_ig).^(-sigma_star) ).*( delta_ig.^(1-sigma_star) ) ;
    deltaX_g_term = sum(deltaX_ig_term,1)';
    deltaX_s = Dsg*deltaX_g_term;
    
    %terms for import tariff revenue
    R_g1 = ( aM_g.^( (1+omega_star)/(1+gammaM*omega_star*eta) ) ).*( varphi_g.^( ( (sigma-eta)*(1+omega_star)  )/( (1-sigma)*(1+gammaM*omega_star*eta) ) ) );
    R_g = R_g1.*varphiR_g;
    R_g(aM_g==0) = 0; %wlog assign zero to products without imports
    R_g(mu_gs_aux<tol) = 0; %wlog assign equal shifters to all products if there are no imports in a sector

    R_s = ( mu_s.^( ( (eta-kappa)*(1+omega_star) )/( (1-eta)*(1+gammaM*omega_star*kappa) ) ) ).*( Dsg*R_g ) ;  
    R_s = DM*R_s;   

    %terms for export tax revenue
    RX_ig_term = a_star_ig.*( (1 + tau_star_ig).^(-sigma_star) ).*( delta_ig.^(1-sigma_star) ).*tau_star_ig ;
    RX_g_term = sum(RX_ig_term,1)';
    RX_s = Dsg*RX_g_term;    
    RX_s = DX*RX_s;

%Compute initial guess of import price index

    PM_s = zeta_s.*( D_s_last.^( gammaM*omega_star/(1+gammaM*omega_star*kappa) ) );
    PM_s(zeta_s==0)=1; %wlog assign price of one to sectors without imports
    %PM_s(m_s0==0)=1; %wlog assign price of one to sectors without imports

%% Solve for prices

% Outter loop: Iterate to solve for sector-level import price
p_s_last = p_s_0;
for j=1:10000
    
    % Inner loop: Iterate to solve for domestic sector-level price (given PM_s)
    p_s=p_s_last;
    for i=1:10000
        %sector-level price index
        P_s = ( AD_s.*( p_s.^(1-kappa) ) + AM_s.*( PM_s.^(1-kappa) ) ).^(1/(1-kappa)) ;
        
        %sector-level input cost
        if nu == 1
            phi_s = exp( alpha_ks'*log( P_s ) );
        else
            phi_s = ( alpha_ks'*( P_s.^(1-nu) ) ).^( 1/(1-nu) );
        end
        
        %sector-region wage   
        w_rs_sectorterm = ( ( alphaL_s.^alphaK_s ).*( p_s ).*( phi_s.^(-alphaI_s) ) ).^(1./(1-alphaI_s)) ;
        w_rs_regionterm = ( ( barL_rs.^(-alphaK_rs) ).*( Z_rs ) ).^(1./(1-alphaI_rs)) ;
        w_rs_regionterm( isnan(w_rs_regionterm) ) = 1; %set wage = 0 for every region-sector without employment
        
        w_rs = w_rs_regionterm.*repmat(w_rs_sectorterm',R,1) ;
        
        %sector-level output value
        Y_s_regionterm = sum( ( Z_rs.*( w_rs.^(-alphaL_rs) ) ).^( 1./alphaK_rs ) ,1)';
        Y_s_sectorterm = ( p_s.*( phi_s.^(-alphaI_s) ) ).^( 1./alphaK_s );
        Y_s = Y_s_sectorterm.*Y_s_regionterm;
        %max(abs( Y_s./sales_s - 1))
        
        %sector-level export value
        X_s = ( deltaX_s ).*( p_s.^(1-sigma_star) );
        %max( abs(X_s./(Dsg*sum(x_ig,1)') - 1))
        
        %Export tax revenue
        RevX = ( RX_s' )*( p_s.^(1-sigma_star) );

        %sector-level final spending shares
        eF_s = ( beta_s.*( P_s.^(1-epsilon) ) )./( beta_s'*( P_s.^(1-epsilon) ) ) ;
        
        %sector-to-sector intermediate spending shares
        eI_ks_sterm = alpha_ks'*( P_s.^(1-nu) );
        eI_ks = diag( P_s.^(1-nu) )*alpha_ks*diag(1./eI_ks_sterm);
        

        %sector-level spending value (final + intermediate)
        tildeE_s = ( eF_s )*(D + RevX + (1-alphaI_s)'*Y_s ) + ( eI_ks*diag(alphaI_s) )*Y_s ;       

        parmR_s = (1 + omega_star)/(1+gammaM*kappa*omega_star);
        R_s_aux = ( ( ( P_s.^(kappa-1) ).*AM_s ).^parmR_s ).*R_s ;
        r_ks = ( repmat( eF_s , 1, S) ).*( repmat( R_s_aux'  , S, 1) );

        E_s_last = ( eye(S)/(eye(S) - r_ks) )*tildeE_s;
        for lcount = 1:1000
            E_s = r_ks*( E_s_last.^parmR_s ) + tildeE_s;

            Edif_s = E_s./E_s_last - 1 ;
            err_Edif_step = max(abs(Edif_s));
            
            if err_Edif_step > tol/100
                E_s_next = E_s_last - adj_outter*(E_s_last - E_s);
                E_s_last = E_s_next;
            else
               break 
            end            
        end    
        
        %sector-level domestic spending value
        ED_s = AD_s.*( (p_s./P_s).^(1-kappa) ).*E_s;
        D_s = E_s.*( P_s.^(kappa-1) );
        
        %sector-level excess supply value
        ESS_s = 1 - (ED_s + X_s)./Y_s ;
        err_ESS_step = max(abs(ESS_s));
        
        %Iterate on domestic sector-level prices
        if err_ESS_step > tol
           p_s_last = p_s;
           p_s = p_s_last - adj_inner*ESS_s;
        else
           break 
        end
    end
    
    %Compute sector-level import price index
    PM_s_next = zeta_s.*( D_s.^( gammaM*omega_star/(1+gammaM*omega_star*kappa) ) );
    PM_s_next(zeta_s==0)=1;
    dif_PM_s = PM_s - PM_s_next ;
    
    err_PM_step = max(abs(dif_PM_s));
    
    %Adjust sector-level import price index
    if err_PM_step > tol
        D_s_next = D_s_last - adj_outter*(D_s_last - D_s);
        
        PM_s = zeta_s.*( D_s_next.^( omega_star/(1+omega_star*kappa) ) );
        PM_s(zeta_s==0)=1;
        D_s_last = D_s_next;
    else
        break 
    end

end

if isnan(err_PM_step) == 1 || isinf(err_PM_step) == 1 || i == 10000
   error('Algorithm did not converge') 
end

%% Compute trade outcomes

%sector-level imports (inclusive of tariffs)
M_s = E_s.*AM_s.*( P_s.^(kappa-1) ).*( PM_s.^(-kappa) );

%product-level imports (inclusive of tariffs)
M_shifter_s = M_s.*( PM_s.^eta );
M_shifter_gs = (M_shifter_s'*Dsg)';

m_g =  ( (   aM_g.*M_shifter_gs                        ).*( varphi_g.^( -eta*(1+gammaM*omega_star*sigma)/(1-sigma) ) ) ).^( 1/(1+gammaM*omega_star*eta) ) ; 
pM_g = ( ( ( aM_g.*M_shifter_gs ).^(gammaM*omega_star) ).*( varphi_g.^(      (1+gammaM*omega_star*sigma)/(1-sigma) ) ) ).^( 1/(1+gammaM*omega_star*eta) ) ; 
pM_g(m_g0 == 0) = 1;

%variety-level imports
M_shifter_g = m_g.*( pM_g.^sigma );
M_shifter_ig = repmat(M_shifter_g', N, 1);

m_ig =      ( (  M_shifter_ig.*a_ig              ).*( z_star_ig.^(-gammaM*sigma) ).*( (1+tau_ig).^(-gammaM*sigma)            ) ).^( 1/(1+gammaM*omega_star*sigma) );
p_star_ig = ( ( (M_shifter_ig.*a_ig).^omega_star ).*( z_star_ig                  ).*( (1+tau_ig).^(-gammaM*sigma*omega_star) ) ).^( 1/(1+gammaM*omega_star*sigma) );

%wlog assign arbitrary values for varieties with zero initial imports
aux = 1./(1+tau_ig);
p_star_ig(m_ig_0==0) = aux(m_ig_0==0); %this is a slow adjustment. I can drop wlog
m_ig(m_ig_0==0) = 0;
pM_ig = (1+tau_ig).*p_star_ig;

%Import values (exclusive tariffs)
VM_ig = m_ig.*p_star_ig;
VM_g = sum(VM_ig,1)';
VM_s = Dsg*VM_g;

%variety-level exports
p_igs = repmat( (p_s'*Dsg) , N, 1);
p_ig = delta_ig.*p_igs;

x_ig = a_star_ig.*( (1+tau_star_ig).*p_ig ).^(-sigma_star);
VX_ig = x_ig.*p_ig;

pX_ig = p_ig.^(gammaX) ; 

%product-level exports
x_g = sum( x_ig,1)';
VX_g = sum(VX_ig,1)';
pX_g = (p_s'*Dsg)';
%pX_g = VX_g./x_g;
%pX_g(x_g==0)=1;

%sector-level exports
VX_s = Dsg*VX_g;
pX_s = (Dsg*pX_g)./sum(Dsg,2);
x_s = VX_s./pX_s;

%sector-level import revenue
RevM_ig = tau_ig.*VM_ig;
RevM_g = sum(RevM_ig, 1)';
RevM_s = Dsg*RevM_g;

clear z_star_ig a_ig w_rs w_rs_regionterm x_ig pX_ig aux pX_igs M_shifter_ig p_star_ig pM_ig p_ig

%% Compute jacobian matrices of outcomes wrt shocks in import costs (zeta) and export demand (delta)

%price-index
eD_s = ED_s./E_s;
psiP_s = 1./( 1 - (1 - eD_s)*( gammaM*omega_star*(kappa-1) )/(1+gammaM*omega_star*kappa)  ) ;
psiPp_ss = diag( psiP_s.*eD_s );
psiPE_ss = diag( ( psiP_s.*(1-eD_s) )*( gammaM*omega_star/(1+gammaM*omega_star*kappa) ) );
psiPz_ss = diag( psiP_s.*(1-eD_s) );

%output
psiYp_ss = diag( 1./(1-alphaI_s) ); 
psiYP_sk = ( diag( alphaI_s./(1-alphaI_s) ) )*( alpha_ks' ) ;

%spending
psiEY_sk = ( repmat( beta_s./E_s , 1, S).*repmat( ( (1-alphaI_s).*Y_s )' , S, 1) ) + ( diag(1./E_s)*alpha_ks*diag( alphaI_s.*Y_s ) );
psiEE_sk = ( repmat( beta_s./E_s , 1, S).*( (1+omega_star)/(1+gammaM*omega_star*kappa) ) ).*repmat( ( ( ( AM_s.*E_s.*( P_s.^(kappa-1) ) ).^( (1+omega_star)/(1+gammaM*omega_star*kappa) ) ).*R_s)' , S, 1);
psiEP_sk = psiEE_sk*(kappa-1);
psiEp_sk = ( repmat( beta_s./E_s , 1, S).*( 1-sigma_star ) ).*(  repmat( ( (p_s.^(1-sigma_star)).*RX_s )' , S, 1) );

%domestic demand
psiDp_ss = diag( (1-kappa)*ED_s./Y_s + (1-sigma_star)*X_s./Y_s );
psiDP_ss = diag( (kappa-1)*ED_s./Y_s );
psiDE_ss = diag( ED_s./Y_s );
psiDd_ss = diag( X_s./Y_s );

%equilibrium price
I = eye(S);
tilde_rhoEeta = ( I/( I - psiEE_sk + (psiEY_sk*psiYP_sk - psiEP_sk )*psiPE_ss ) );
tilde_rhoEp = tilde_rhoEeta*( psiEp_sk + psiEP_sk*psiPp_ss + psiEY_sk*(psiYp_ss - psiYP_sk*psiPp_ss) );
tilde_rhoEz = tilde_rhoEeta*( (psiEP_sk - psiEY_sk*psiYP_sk)*psiPz_ss );

rhop = I/( psiYp_ss - psiDp_ss - (psiYP_sk + psiDP_ss)*psiPp_ss - ( (psiYP_sk + psiDP_ss)*psiPE_ss + psiDE_ss)*tilde_rhoEp  );
rho_pzeta = rhop*( (psiYP_sk+psiDP_ss)*(psiPz_ss + psiPE_ss*tilde_rhoEz) + psiDE_ss*tilde_rhoEz );
rho_pdelta = rhop*psiDd_ss;
rho_peta = rhop*( (psiYP_sk+psiDP_ss)*psiPE_ss + psiDE_ss )*tilde_rhoEeta;

%% Compute jacobian of import cost (zeta) and export demand (delta) wrt tariffs 
%sectors in rows, varieties ig in columns (ordered first by country and then by product)

%Jacobian of zeta
rho_zeta_term_g = ( mu_g./( ( mu_s'*Dsg )' ) )./varphi_g ;
rho_zeta_term_ig = ( 1/(1+gammaM*omega_star*kappa) )*( varphi_ig*spdiags( rho_zeta_term_g , 0 , G, G) ) ;
rho_zeta_term_ig(:,mu_gs_aux' < tol) = 0; %assign zero effect of varieties in sectors without imports
%rho_zeta_term_g(m_ig_0==0) = 0;

rho_zeta_tau_vec_ig     = reshape( (rho_zeta_term_ig.*( gammaM./(1+tau_ig)             ) )', N*G, 1); 
rho_zeta_tau_s_ig       = repmat( Dsg, 1, N).*repmat( rho_zeta_tau_vec_ig', S, 1);
rho_zeta_tau_s_ig       = sparse( rho_zeta_tau_s_ig   )*rhoM;

rho_zeta_zstar_vec_ig   = reshape( (rho_zeta_term_ig.*( gammaM                         ) )', N*G, 1); 
rho_zeta_zstar_s_ig     = repmat( Dsg, 1, N).*repmat( rho_zeta_zstar_vec_ig'  , S, 1);
rho_zeta_zstar_s_ig     = sparse( rho_zeta_zstar_s_ig )     ;

rho_zeta_a_vec_ig       = reshape( (rho_zeta_term_ig.*( (1+gammaM*omega_star)/(1-sigma) ))', N*G, 1); 
rho_zeta_a_s_ig         = repmat( Dsg, 1, N).*repmat( rho_zeta_a_vec_ig'  , S, 1);
rho_zeta_a_s_ig         = sparse( rho_zeta_a_s_ig     )     ;


clear rho_zeta_term_g rho_zeta_term_ig rho_zeta_tau_vec_ig rho_zeta_zstar_vec_ig rho_zeta_a_vec_ig

%Jacobian of delta
rho_delta_term_s = 1./( deltaX_s )';
rho_delta_term_s( isinf(rho_delta_term_s) ) = 0; %assign zero effect on delta for varieties in sectors without exports
rho_delta_term_g = ( rho_delta_term_s*Dsg )';


rho_delta_term_ig = - (rhoX*sigma_star)*( deltaX_ig_term*spdiags( rho_delta_term_g , 0, G, G  ) )./(1+tau_star_ig);
rho_delta_vec_ig = reshape( rho_delta_term_ig', N*G, 1); 
rho_delta_taustar_s_ig = repmat( Dsg, 1, N).*repmat( rho_delta_vec_ig', S, 1);
rho_delta_taustar_s_ig = sparse( rho_delta_taustar_s_ig );

rho_delta_term_ig =                     ( deltaX_ig_term*spdiags( rho_delta_term_g , 0, G, G  ) ) ;
rho_delta_vec_ig = reshape( rho_delta_term_ig', N*G, 1); 
rho_delta_astar_s_ig = repmat( Dsg, 1, N).*repmat( rho_delta_vec_ig', S, 1);
rho_delta_astar_s_ig = sparse( rho_delta_astar_s_ig );

rho_delta_term_ig = (1-sigma_star)     *( deltaX_ig_term*spdiags( rho_delta_term_g , 0, G, G  ) ) ;
rho_delta_vec_ig = reshape( rho_delta_term_ig', N*G, 1); 
rho_delta_delta_s_ig = repmat( Dsg, 1, N).*repmat( rho_delta_vec_ig', S, 1);
rho_delta_delta_s_ig = sparse( rho_delta_delta_s_ig );

clear rho_delta_term_s  rho_delta_term_g rho_delta_term_ig rho_delta_vec_ig deltaX_ig_term 

%Jacobian of eta
%terms associated with import tariff revenue
etaM_term1_s = ( ( AM_s.*E_s.*( P_s.^(kappa-1) ) ).^( (1+omega_star)/(1+gammaM*omega_star*kappa) ) ).*( R_s./mu_s )*( ( (eta-kappa)*(1+omega_star) )/( (1+gammaM*omega_star*kappa)*(1+gammaM*omega_star*eta) ) );
etaM_term1_s(  Dsg*m_g0<tol ) = 0;
etaM_term1_g = (etaM_term1_s'*Dsg)'.*(mu_g./varphi_g);
etaM_term1_ig = ( varphi_ig )*spdiags( etaM_term1_g,  0, G, G  );

etaM_term2_s = ( ( AM_s.*E_s.*( P_s.^(kappa-1) ) ).^( (1+omega_star)/(1+gammaM*omega_star*kappa) ) ).*( mu_s.^( ( (eta-kappa)*(1+omega_star) )/( (1-eta)*(1+gammaM*omega_star*kappa) ) ) ) ;
etaM_term2_s(  Dsg*m_g0<tol ) = 0;
etaM_term2_g = ( (etaM_term2_s'*Dsg )'.*R_g1.*varphiR_g./varphi_g )*( ( (sigma-eta)*(1+omega_star) )/( (1+gammaM*omega_star*sigma)*(1+gammaM*omega_star*eta) ) );
etaM_term2_ig = ( varphi_ig )*spdiags( etaM_term2_g,  0, G, G  );

etaM_term3_s = etaM_term2_s ;
etaM_term3_g = (etaM_term3_s'*Dsg )'.*R_g1;
etaM_term3_ig = tilde_varphiR_ig*spdiags( etaM_term3_g ,  0, G, G  );


rho_etaM_term1_ig = ( etaM_term1_ig + etaM_term2_ig ).*( gammaM./(1+tau_ig) ) ;
%rho_etaM_term3_ig = etaM_term3_ig.*( 1 - ( (1+omega_star)*sigma*gammaM/(1+omega_star*sigma*gammaM) ).*( tau_ig./(1+tau_ig) )  ) ;
rho_etaM_term3_ig = etaM_term3_ig.*( gammaM - ( (1+omega_star)*sigma*gammaM/(1+omega_star*sigma*gammaM) ).*( tau_ig./(1+tau_ig) )  ) ;
rho_etaM_ig = DM*( rho_etaM_term1_ig + rho_etaM_term3_ig );
rho_etaM_vec_ig = reshape( rho_etaM_ig', N*G, 1); 
rho_eta_tau_s_ig = repmat( beta_s./E_s , 1, N*G).*repmat( rho_etaM_vec_ig', S, 1);
rho_eta_tau_s_ig = sparse( rho_eta_tau_s_ig )*rhoM ;

rho_etaM_term1_ig = ( etaM_term1_ig + etaM_term2_ig ).*gammaM ;
rho_etaM_term3_ig = etaM_term3_ig.*( (1 - sigma*gammaM)/(1+omega_star*sigma*gammaM) ).*( tau_ig ) ;
rho_etaM_ig = DM*( rho_etaM_term1_ig + rho_etaM_term3_ig );
rho_etaM_vec_ig = reshape( rho_etaM_ig', N*G, 1); 
rho_eta_zstar_s_ig = repmat( beta_s./E_s , 1, N*G).*repmat( rho_etaM_vec_ig', S, 1);
rho_eta_zstar_s_ig = sparse( rho_eta_zstar_s_ig ) ;

rho_etaM_term1_ig = ( etaM_term1_ig + etaM_term2_ig ).*( (1+omega_star*gammaM)/(1-sigma) );
rho_etaM_term3_ig = etaM_term3_ig.*( (1 + omega_star)/(1+omega_star*sigma*gammaM) ).*( tau_ig ) ;
rho_etaM_ig = DM*( rho_etaM_term1_ig + rho_etaM_term3_ig );
rho_etaM_vec_ig = reshape( rho_etaM_ig', N*G, 1); 
rho_eta_a_s_ig = repmat( beta_s./E_s , 1, N*G).*repmat( rho_etaM_vec_ig', S, 1);
rho_eta_a_s_ig = sparse( rho_eta_a_s_ig ) ;


clear etaM_term1_ig etaM_term2_ig etaM_term3_ig etaM_term1_g etaM_term2_g etaM_term3_g etaM_term1_s etaM_term2_s etaM_term3_s rho_etaM_term1_ig rho_etaM_term3_ig rho_etaM_ig rho_etaM_vec_ig

%terms associated with export tax revenue
etaX_term1_s = - sigma_star*( p_s.^(1-sigma_star) );
etaX_term1_ig =  RX_ig_term./(1+tau_star_ig) ;
etaX_term1_ig = etaX_term1_ig*spdiags( (etaX_term1_s'*Dsg )' ,  0, G, G  );

etaX_term2_s =  p_s.^(1-sigma_star) ;
etaX_term2_ig = a_star_ig.*( (1 + tau_star_ig).^(-sigma_star) ).*( delta_ig.^(1-sigma_star) ) ;
etaX_term2_ig = etaX_term2_ig*spdiags( (etaX_term2_s'*Dsg )' ,  0, G, G  );

rho_etaX_ig = DX*( etaX_term1_ig + etaX_term2_ig );
rho_etaX_vec_ig = reshape( rho_etaX_ig', N*G, 1); 
rho_eta_taustar_s_ig = repmat( beta_s./E_s , 1, N*G).*repmat( rho_etaX_vec_ig', S, 1);
rho_eta_taustar_s_ig = sparse( rho_eta_taustar_s_ig );

if DX == 1
    adj_etaX = spdiags( tau_star_ig./( 1 - sigma_star*tau_star_ig./(1+tau_star_ig)  ) , 0, N*G, N*G);
else
    adj_etaX = 0;
end

rho_eta_astar_s_ig = rho_eta_taustar_s_ig*adj_etaX;
rho_eta_delta_s_ig = rho_eta_astar_s_ig*(1-sigma_star);
rho_eta_taustar_s_ig = rho_eta_taustar_s_ig*rhoX;

clear etaX_term1_ig etaX_term2_ig etaX_term2_s rho_etaX_ig rho_etaX_vec_ig a_star_ig delta_ig

%Select sample columns of jacobian matrices

rho_zeta_tau_s_ig = rho_zeta_tau_s_ig(:,ind_Mshifts_vec');  
rho_zeta_a_s_ig = rho_zeta_a_s_ig(:,ind_Mshifts_vec');  
rho_zeta_zstar_s_ig = rho_zeta_zstar_s_ig(:,ind_Mshifts_vec');  

rho_eta_tau_s_ig = rho_eta_tau_s_ig(:,ind_Mshifts_vec'); 
rho_eta_a_s_ig = rho_eta_a_s_ig(:,ind_Mshifts_vec'); 
rho_eta_zstar_s_ig = rho_eta_zstar_s_ig(:,ind_Mshifts_vec'); 

rho_delta_taustar_s_ig = rho_delta_taustar_s_ig(:,ind_Xshifts_vec');  
rho_delta_astar_s_ig = rho_delta_astar_s_ig(:,ind_Xshifts_vec');  
rho_delta_delta_s_ig = rho_delta_delta_s_ig(:,ind_Xshifts_vec');  

rho_eta_taustar_s_ig = rho_eta_taustar_s_ig(:,ind_Xshifts_vec'); 
rho_eta_astar_s_ig = rho_eta_astar_s_ig(:,ind_Xshifts_vec'); 
rho_eta_delta_s_ig = rho_eta_delta_s_ig(:,ind_Xshifts_vec'); 


%% Jacobian of sector-level outcomes wrt fundamental shocks

% Jacobian of prices
rho_p_tau_s_ig   =     sparse( rho_pzeta*rho_zeta_tau_s_ig   + rho_peta*rho_eta_tau_s_ig   );
rho_p_a_s_ig     =     sparse( rho_pzeta*rho_zeta_a_s_ig     + rho_peta*rho_eta_a_s_ig     );
rho_p_zstar_s_ig =     sparse( rho_pzeta*rho_zeta_zstar_s_ig + rho_peta*rho_eta_zstar_s_ig );

rho_p_taustar_s_ig = sparse( rho_pdelta*rho_delta_taustar_s_ig + rho_peta*rho_eta_taustar_s_ig );
rho_p_astar_s_ig = sparse( rho_pdelta*rho_delta_astar_s_ig     + rho_peta*rho_eta_astar_s_ig );
rho_p_delta_s_ig = sparse( rho_pdelta*rho_delta_delta_s_ig     + rho_peta*rho_eta_delta_s_ig );

%Jacobian of spending 
rho_Ezeta = tilde_rhoEp*rho_pzeta + tilde_rhoEz ; 
rho_Edelta = tilde_rhoEp*rho_pdelta ;
rho_Eeta =   tilde_rhoEp*rho_peta  + tilde_rhoEeta;

rho_E_tau_s_ig   =     sparse( rho_Ezeta*rho_zeta_tau_s_ig   + rho_Eeta*rho_eta_tau_s_ig   );
rho_E_a_s_ig     =     sparse( rho_Ezeta*rho_zeta_a_s_ig     + rho_Eeta*rho_eta_a_s_ig     );
rho_E_zstar_s_ig =     sparse( rho_Ezeta*rho_zeta_zstar_s_ig + rho_Eeta*rho_eta_zstar_s_ig );

rho_E_taustar_s_ig = sparse( rho_Edelta*rho_delta_taustar_s_ig   + rho_Eeta*rho_eta_taustar_s_ig );
rho_E_astar_s_ig   = sparse( rho_Edelta*rho_delta_astar_s_ig     + rho_Eeta*rho_eta_astar_s_ig   );
rho_E_delta_s_ig   = sparse( rho_Edelta*rho_delta_delta_s_ig     + rho_Eeta*rho_eta_delta_s_ig   );

%Jacobian of price index
rho_Pzeta =  psiPp_ss*rho_pzeta  + psiPE_ss*rho_Ezeta + psiPz_ss ; 
rho_Pdelta = psiPp_ss*rho_pdelta + psiPE_ss*rho_Edelta;
rho_Peta =   psiPp_ss*rho_peta   + psiPE_ss*rho_Eeta ;

rho_P_tau_s_ig   =     sparse( rho_Pzeta*rho_zeta_tau_s_ig   + rho_Peta*rho_eta_tau_s_ig   );
rho_P_a_s_ig     =     sparse( rho_Pzeta*rho_zeta_a_s_ig     + rho_Peta*rho_eta_a_s_ig     );
rho_P_zstar_s_ig =     sparse( rho_Pzeta*rho_zeta_zstar_s_ig + rho_Peta*rho_eta_zstar_s_ig );

rho_P_taustar_s_ig = sparse( rho_Pdelta*rho_delta_taustar_s_ig   + rho_Peta*rho_eta_taustar_s_ig );
rho_P_astar_s_ig   = sparse( rho_Pdelta*rho_delta_astar_s_ig     + rho_Peta*rho_eta_astar_s_ig   );
rho_P_delta_s_ig   = sparse( rho_Pdelta*rho_delta_delta_s_ig     + rho_Peta*rho_eta_delta_s_ig   );

%Jacobian of output 
% rho_Yzeta =  psiYp_ss*rho_pzeta  - psiYP_sk*rho_Pzeta; 
% rho_Ydelta = psiYp_ss*rho_pdelta - psiYP_sk*rho_Pdelta;
% rho_Yeta =   psiYp_ss*rho_peta   - psiYP_sk*rho_Peta ;
% 
% rho_Y_tau_s_ig   =     sparse( rho_Yzeta*rho_zeta_tau_s_ig   + rho_Yeta*rho_eta_tau_s_ig );
% rho_Y_a_s_ig     =     sparse( rho_Yzeta*rho_zeta_a_s_ig     + rho_Yeta*rho_eta_a_s_ig     );
% rho_Y_zstar_s_ig =     sparse( rho_Yzeta*rho_zeta_zstar_s_ig + rho_Yeta*rho_eta_zstar_s_ig );
% 
% rho_Y_taustar_s_ig = sparse( rho_Ydelta*rho_delta_taustar_s_ig   + rho_Yeta*rho_eta_taustar_s_ig );
% rho_Y_astar_s_ig   = sparse( rho_Ydelta*rho_delta_astar_s_ig     + rho_Yeta*rho_eta_astar_s_ig   );
% rho_Y_delta_s_ig   = sparse( rho_Ydelta*rho_delta_delta_s_ig     + rho_Yeta*rho_eta_delta_s_ig   );

%Jacobian of import price index 
rho_PMzeta  =  ( omega_star*gammaM/(1+omega_star*kappa*gammaM) )*( rho_Ezeta  + (kappa-1)*rho_Pzeta ) + I ; 
rho_PMdelta =  ( omega_star*gammaM/(1+omega_star*kappa*gammaM) )*( rho_Edelta + (kappa-1)*rho_Pdelta);
rho_PMeta   =  ( omega_star*gammaM/(1+omega_star*kappa*gammaM) )*( rho_Eeta   + (kappa-1)*rho_Peta  );

rho_PM_tau_s_ig   =     sparse( rho_PMzeta*rho_zeta_tau_s_ig   + rho_PMeta*rho_eta_tau_s_ig );
rho_PM_a_s_ig     =     sparse( rho_PMzeta*rho_zeta_a_s_ig     + rho_PMeta*rho_eta_a_s_ig     );
rho_PM_zstar_s_ig =     sparse( rho_PMzeta*rho_zeta_zstar_s_ig + rho_PMeta*rho_eta_zstar_s_ig );

rho_PM_taustar_s_ig = sparse( rho_PMdelta*rho_delta_taustar_s_ig   + rho_PMeta*rho_eta_taustar_s_ig );
rho_PM_astar_s_ig   = sparse( rho_PMdelta*rho_delta_astar_s_ig     + rho_PMeta*rho_eta_astar_s_ig   );
rho_PM_delta_s_ig   = sparse( rho_PMdelta*rho_delta_delta_s_ig     + rho_PMeta*rho_eta_delta_s_ig   );

%Jacobian of export value wrt tariff
% rho_Xzeta  =  (1-sigma_star)*rho_pzeta; 
% rho_Xdelta =  (1-sigma_star)*rho_pdelta + I ;
% rho_Xeta   =  (1-sigma_star)*rho_peta;
% 
% rho_vx_tau_s_ig   =     sparse( rho_Xzeta*rho_zeta_tau_s_ig   + rho_Xeta*rho_eta_tau_s_ig );
% rho_vx_a_s_ig     =     sparse( rho_Xzeta*rho_zeta_a_s_ig     + rho_Xeta*rho_eta_a_s_ig     );
% rho_vx_zstar_s_ig =     sparse( rho_Xzeta*rho_zeta_zstar_s_ig + rho_Xeta*rho_eta_zstar_s_ig );
% 
% rho_vx_taustar_s_ig = sparse( rho_Xdelta*rho_delta_taustar_s_ig   + rho_Xeta*rho_eta_taustar_s_ig );
% rho_vx_astar_s_ig   = sparse( rho_Xdelta*rho_delta_astar_s_ig     + rho_Xeta*rho_eta_astar_s_ig   );
% rho_vx_delta_s_ig   = sparse( rho_Xdelta*rho_delta_delta_s_ig     + rho_Xeta*rho_eta_delta_s_ig   );


%Jacobian of imports (inclusinve of tariffs) wrt tariffs 
% rho_VMzeta  =  (1-kappa)*(rho_PMzeta  - rho_Pzeta ) + rho_Ezeta  ; 
% rho_VMdelta =  (1-kappa)*(rho_PMdelta - rho_Pdelta) + rho_Edelta ;
% rho_VMeta   =  (1-kappa)*(rho_PMeta   - rho_Peta  ) + rho_Eeta   ;
% 
% rho_M_tau_s_ig   =     sparse( rho_VMzeta*rho_zeta_tau_s_ig   + rho_VMeta*rho_eta_tau_s_ig );
% rho_M_a_s_ig     =     sparse( rho_VMzeta*rho_zeta_a_s_ig     + rho_VMeta*rho_eta_a_s_ig     );
% rho_M_zstar_s_ig =     sparse( rho_VMzeta*rho_zeta_zstar_s_ig + rho_VMeta*rho_eta_zstar_s_ig );
% 
% rho_M_taustar_s_ig = sparse( rho_VMdelta*rho_delta_taustar_s_ig   + rho_VMeta*rho_eta_taustar_s_ig );
% rho_M_astar_s_ig   = sparse( rho_VMdelta*rho_delta_astar_s_ig     + rho_VMeta*rho_eta_astar_s_ig   );
% rho_M_delta_s_ig   = sparse( rho_VMdelta*rho_delta_delta_s_ig     + rho_VMeta*rho_eta_delta_s_ig   );

clear etaM_term1_ig etaM_term2_ig etaM_term3_ig rho_VM_ig rho_VM_vec_ig
clear rho_VMzeta rho_VMdelta rho_VMeta
clear rho_Xzeta  rho_Xdelta rho_Xeta
clear rho_PMzeta  rho_PMdelta rho_PMeta
clear rho_Yzeta  rho_Ydelta rho_Yeta
clear rho_Pzeta  rho_Pdelta rho_Peta
clear rho_Ezeta  rho_Edelta rho_Eeta

clear rho_eta_taustar_s_ig rho_delta_taustar_s_ig rho_eta_astar_s_ig rho_delta_astar_s_ig rho_eta_delta_s_ig rho_delta_delta_s_ig 
clear rho_eta_tau_s_ig rho_zeta_tau_s_ig rho_eta_a_s_ig rho_zeta_a_s_ig rho_eta_zstar_s_ig rho_zeta_zstar_s_ig
clear rho_etaM_vec_ig m_ig tau

%% Jacobian for variety-level outcomes

%Jacobian of import quantity
    %term 1 - variety
    ind_ig = m_ig_0 > 0;
    tau_ig_vec = reshape( ( ind_ig./(1 + tau_ig) )', N*G, 1); 
    rho_m_ig_tau_ig_term1 = -( sigma*gammaM/(1+omega_star*sigma*gammaM) )*spdiags(tau_ig_vec, 0, N*G, N*G);
    rho_m_ig_tau_ig_term1   = rhoM*rho_m_ig_tau_ig_term1(ind_Msample_ig_vec,ind_Mshifts_vec');

    rho_m_ig_a_ig_term1 = (     1/(1+omega_star*sigma*gammaM) )*speye(N*G); %spdiags(ind_ig_vec, 0, N*G, N*G);
    rho_m_ig_a_ig_term1    =       rho_m_ig_a_ig_term1(ind_Msample_ig_vec,ind_Mshifts_vec');
    
    rho_m_ig_zstar_ig_term1 = -( sigma*gammaM/(1+omega_star*sigma*gammaM) )*speye(N*G); %*spdiags(ind_ig_vec, 0, N*G, N*G);
    rho_m_ig_zstar_ig_term1 =      rho_m_ig_zstar_ig_term1(ind_Msample_ig_vec,ind_Mshifts_vec');

    %term 2 - product
    den_rho_m_ig_tau_ig_term2 = 1./varphi_g;
    den_rho_m_ig_tau_ig_term2(m_g0<tol) = 0;
    
    aux_rho_m_ig_tau_ig_term2   = ( (sigma-eta)/( (1+omega_star*eta*gammaM)*(1+omega_star*sigma*gammaM) ) )*( ind_ig.*varphi_ig.*gammaM./(1 + tau_ig)            )*spdiags( den_rho_m_ig_tau_ig_term2 , 0, G, G);
    aux_rho_m_ig_a_ig_term2     = ( (sigma-eta)/( (1+omega_star*eta*gammaM)*(1+omega_star*sigma*gammaM) ) )*( ind_ig.*varphi_ig*(1+ gammaM*omega_star)/(1-sigma) )*spdiags( den_rho_m_ig_tau_ig_term2 , 0, G, G);
    aux_rho_m_ig_zstar_ig_term2 = ( (sigma-eta)/( (1+omega_star*eta*gammaM)*(1+omega_star*sigma*gammaM) ) )*( ind_ig.*varphi_ig.*gammaM                          )*spdiags( den_rho_m_ig_tau_ig_term2 , 0, G, G);

    aux_rho_m_ig_tau_ig_term2   = sparse( reshape(aux_rho_m_ig_tau_ig_term2'  , N*G, 1) );
    aux_rho_m_ig_a_ig_term2     = sparse( reshape(aux_rho_m_ig_a_ig_term2'    , N*G, 1) );
    aux_rho_m_ig_zstar_ig_term2 = sparse( reshape(aux_rho_m_ig_zstar_ig_term2', N*G, 1) );

    aux_rho_m_ig_tau_ig_term2   = rhoM*aux_rho_m_ig_tau_ig_term2(ind_Mshifts_vec)' ;
    aux_rho_m_ig_a_ig_term2     =      aux_rho_m_ig_a_ig_term2(ind_Mshifts_vec)' ;    
    aux_rho_m_ig_zstar_ig_term2 =      aux_rho_m_ig_zstar_ig_term2(ind_Mshifts_vec)' ;


    good_ind = repmat(spdiags(ones(G,1), 0, G, G), 1, N);
    good_ind = good_ind(:, ind_Mshifts_vec');
    
    rho_m_ig_tau_ig_term2   = [];
    rho_m_ig_a_ig_term2     = [];
    rho_m_ig_zstar_ig_term2 = [];
    Nstep = floor( G/2000 ) + 1;
    for g = 1:Nstep
        g1 = (g-1)*2000+1;
        g2 = g*2000;
        if g2 > G
            g2 = G;
        end
        aux2_rho_m_ig_tau_ig_term2   = good_ind(g1:g2,:).*repmat(aux_rho_m_ig_tau_ig_term2, g2-g1+1, 1 );
        aux2_rho_m_ig_a_ig_term2     = good_ind(g1:g2,:).*repmat(aux_rho_m_ig_a_ig_term2, g2-g1+1, 1 );
        aux2_rho_m_ig_zstar_ig_term2 = good_ind(g1:g2,:).*repmat(aux_rho_m_ig_zstar_ig_term2, g2-g1+1, 1 );

        rho_m_ig_tau_ig_term2   = [rho_m_ig_tau_ig_term2; aux2_rho_m_ig_tau_ig_term2];
        rho_m_ig_a_ig_term2     = [rho_m_ig_a_ig_term2; aux2_rho_m_ig_a_ig_term2];
        rho_m_ig_zstar_ig_term2 = [rho_m_ig_zstar_ig_term2; aux2_rho_m_ig_zstar_ig_term2];

        rho_m_ig_tau_ig_term2   = sparse(rho_m_ig_tau_ig_term2);
        rho_m_ig_a_ig_term2     = sparse(rho_m_ig_a_ig_term2);
        rho_m_ig_zstar_ig_term2 = sparse(rho_m_ig_zstar_ig_term2);
    end    
    rho_m_ig_tau_ig_term2 = repmat(rho_m_ig_tau_ig_term2, N, 1);
    rho_m_ig_a_ig_term2 = repmat(rho_m_ig_a_ig_term2, N, 1);
    rho_m_ig_zstar_ig_term2 = repmat(rho_m_ig_zstar_ig_term2, N, 1);

    rho_m_ig_tau_ig_term2 = rho_m_ig_tau_ig_term2(ind_Msample_ig_vec,:);
    rho_m_ig_a_ig_term2 = rho_m_ig_a_ig_term2(ind_Msample_ig_vec,:);
    rho_m_ig_zstar_ig_term2 = rho_m_ig_zstar_ig_term2(ind_Msample_ig_vec,:);

    %term 3 - sector   
    rho_m_s_tau_ig_term3   =     ( 1/(1+omega_star*eta*gammaM) )*(  rho_E_tau_s_ig   + (kappa-1)*rho_P_tau_s_ig   + (eta-kappa)*rho_PM_tau_s_ig    );
    rho_m_s_a_ig_term3     =     ( 1/(1+omega_star*eta*gammaM) )*(  rho_E_a_s_ig     + (kappa-1)*rho_P_a_s_ig     + (eta-kappa)*rho_PM_a_s_ig      );
    rho_m_s_zstar_ig_term3 =     ( 1/(1+omega_star*eta*gammaM) )*(  rho_E_zstar_s_ig + (kappa-1)*rho_P_zstar_s_ig + (eta-kappa)*rho_PM_zstar_s_ig  );

    dum_ig_sec = repmat( Dsg, 1, N)';
    dum_ig_sec = dum_ig_sec(ind_Msample_ig_vec, :);

    rho_m_ig_tau_ig_term3   = dum_ig_sec*rho_m_s_tau_ig_term3;
    rho_m_ig_a_ig_term3     = dum_ig_sec*rho_m_s_a_ig_term3;
    rho_m_ig_zstar_ig_term3 = dum_ig_sec*rho_m_s_zstar_ig_term3;

    %rho_qm_ig_tau_ig    = (rho_m_ig_tau_ig_term1  +rho_m_ig_tau_ig_term2   +rho_m_ig_tau_ig_term3);
    rho_qm_ig_tau_ig    = (rho_m_ig_tau_ig_term1  +rho_m_ig_tau_ig_term2   +rho_m_ig_tau_ig_term3)*gammaQ;
    rho_qm_ig_a_ig      = rho_m_ig_a_ig_term1    +rho_m_ig_a_ig_term2     +rho_m_ig_a_ig_term3;
    rho_qm_ig_zstar_ig  = rho_m_ig_zstar_ig_term1+rho_m_ig_zstar_ig_term2 +rho_m_ig_zstar_ig_term3;

    %term 4 - sector jacobian wrt taustar    
    rho_m_s_taustar_ig =     ( 1/(1+omega_star*kappa*gammaM) )*( rho_E_taustar_s_ig     + (kappa-1)*rho_P_taustar_s_ig     + (eta-kappa)*rho_PM_taustar_s_ig  );
    rho_m_s_astar_ig   =     ( 1/(1+omega_star*kappa*gammaM) )*( rho_E_astar_s_ig       + (kappa-1)*rho_P_astar_s_ig       + (eta-kappa)*rho_PM_astar_s_ig    );
    rho_m_s_delta_ig   =     ( 1/(1+omega_star*kappa*gammaM) )*( rho_E_delta_s_ig       + (kappa-1)*rho_P_delta_s_ig       + (eta-kappa)*rho_PM_delta_s_ig    );

    %rho_qm_ig_taustar_ig = ( dum_ig_sec*rho_m_s_taustar_ig );
    rho_qm_ig_taustar_ig = ( dum_ig_sec*rho_m_s_taustar_ig )*gammaQ;
    rho_qm_ig_astar_ig = dum_ig_sec*rho_m_s_astar_ig;
    rho_qm_ig_delta_ig = dum_ig_sec*rho_m_s_delta_ig;

    clear dum_ig_sec rho_m_ig_tau_ig_term1   rho_m_ig_tau_ig_term2   rho_m_ig_tau_ig_term3   aux2_rho_m_ig_tau_ig_term2 rho_m_s_taustar_ig;
    clear dum_ig_sec rho_m_ig_a_ig_term1     rho_m_ig_a_ig_term2     rho_m_ig_a_ig_term3     aux2_rho_m_ig_a_ig_term2  rho_m_s_astar_ig;
    clear dum_ig_sec rho_m_ig_zstar_ig_term1 rho_m_ig_zstar_ig_term2 rho_m_ig_zstar_ig_term3 aux2_rho_m_ig_zstar_ig_term2  rho_m_s_delta_ig;

% Jacobian for variety-level import price (pre-tariff)
    rho_pstar_ig_tau_ig     = omega_star*rho_qm_ig_tau_ig    ;
    rho_pstar_ig_a_ig       = omega_star*rho_qm_ig_a_ig    ;

    rho_pstar_ig_zstar_ig_term1 = speye(N*G);
    rho_pstar_ig_zstar_ig_term1 = rho_pstar_ig_zstar_ig_term1(ind_Msample_ig_vec,ind_Mshifts_vec');

    rho_pstar_ig_zstar_ig   = omega_star*rho_qm_ig_zstar_ig  + rho_pstar_ig_zstar_ig_term1;

    rho_pstar_ig_taustar_ig = omega_star*rho_qm_ig_taustar_ig;
    rho_pstar_ig_astar_ig   = omega_star*rho_qm_ig_astar_ig;
    rho_pstar_ig_delta_ig   = omega_star*rho_qm_ig_delta_ig;

% Jacobian for variety-level import price (post-tariff)
    tau_ig_vec = reshape( ( 1./(1 + tau_ig) )', N*G, 1); 
    rho_pm_ig_tau_ig_term1 = spdiags(tau_ig_vec, 0, N*G, N*G);
    rho_pm_ig_tau_ig_term1 = rhoM*rho_pm_ig_tau_ig_term1(ind_Msample_ig_vec,ind_Mshifts_vec');
    %rho_pm_ig_tau_ig_term1 = rhoM*rho_pm_ig_tau_ig_term1(ind_Msample_ig_vec,ind_Mshifts_vec')*gammaQ;

    rho_pm_ig_tau_ig   =  rho_pstar_ig_tau_ig      + rho_pm_ig_tau_ig_term1;
    rho_pm_ig_a_ig     =  rho_pstar_ig_a_ig      ;
    rho_pm_ig_zstar_ig =  rho_pstar_ig_zstar_ig  ;
    
    rho_pm_ig_taustar_ig = rho_pstar_ig_taustar_ig;
    rho_pm_ig_astar_ig   = rho_pstar_ig_astar_ig;
    rho_pm_ig_delta_ig   = rho_pstar_ig_delta_ig;

% Jacobian for variety-level export price        
    dum_ig_sec = repmat( Dsg, 1, N)';
    dum_ig_sec = dum_ig_sec(ind_Xsample_ig_vec, :);

    
    rho_px_ig_delta_ig_term1 = speye(N*G);
    rho_px_ig_delta_ig_term1 = rho_px_ig_delta_ig_term1(ind_Xsample_ig_vec,ind_Xshifts_vec');   

    rho_px_ig_tau_ig   = gammaX*dum_ig_sec*rho_p_tau_s_ig;
    rho_px_ig_a_ig     = gammaX*dum_ig_sec*rho_p_a_s_ig;
    rho_px_ig_zstar_ig = gammaX*dum_ig_sec*rho_p_zstar_s_ig;

    rho_px_ig_taustar_ig = gammaX*dum_ig_sec*rho_p_taustar_s_ig;
    rho_px_ig_astar_ig   = gammaX*dum_ig_sec*rho_p_astar_s_ig;
    rho_px_ig_delta_ig   = gammaX*dum_ig_sec*rho_p_delta_s_ig   +  rho_px_ig_delta_ig_term1 ;

%Jacobian for variety-level revenue term
    rm_tau_ig_vec = reshape( ( tau_ig./(1 + tau_ig) )', N*G, 1); 
    rho_rm_ig_tau_ig_term1 = spdiags(rm_tau_ig_vec, 0, N*G, N*G);
    rho_rm_ig_tau_ig_term1 = rho_rm_ig_tau_ig_term1(ind_Msample_ig_vec,ind_Msample_ig_vec');

    rho_rm_ig_tau_ig   = rho_pm_ig_tau_ig_term1 + rho_rm_ig_tau_ig_term1*(rho_pstar_ig_tau_ig   + rho_qm_ig_tau_ig  );
    rho_rm_ig_a_ig     =                        + rho_rm_ig_tau_ig_term1*(rho_pstar_ig_a_ig     + rho_qm_ig_a_ig    );
    rho_rm_ig_zstar_ig =                        + rho_rm_ig_tau_ig_term1*(rho_pstar_ig_zstar_ig + rho_qm_ig_zstar_ig);

    rho_rm_ig_taustar_ig =                      + rho_rm_ig_tau_ig_term1*(rho_pstar_ig_taustar_ig + rho_qm_ig_taustar_ig);
    rho_rm_ig_astar_ig =                        + rho_rm_ig_tau_ig_term1*(rho_pstar_ig_astar_ig   + rho_qm_ig_astar_ig);
    rho_rm_ig_delta_ig =                        + rho_rm_ig_tau_ig_term1*(rho_pstar_ig_delta_ig   + rho_qm_ig_delta_ig);
   
end

